A ideia do notebook é a partir do dos contratos futuros de DI construir a estrutura a termo de juros da economia brasileira. Por fim, vamos aplicar uma Análise de Componentes Principais e encontrar os fatores de risco latentes na curva de juros. Se você já ouviu falar em Nelson-Siegel, você já deve imaginar quais são.

Na análise, serão usados os ajustes dos contratos negociados na B3 — e tem bastante coisa, de boi-gordo a câmbio, de cupom cambial a café. Os dados podem ser obtidos facilmente por meio do pacote rb3. Vou utilizar também o pacote fixedincome. Ambos foram desenvolvidos pelo Wilson Freitas e facilitam bastante a vida.

Obtendo e transformando os dados da B3

library(tidyverse)
library(rb3) # Dados da B3 
library(fixedincome) # Manipulação de dados de renda fixa e juros
library(bizdays) # Calendários
library(plotly)
library(xts)
library(reactable)
library(psych)

paleta_cores = c('#b8acd1', '#e2d5f1', '#c9e8fd', '#526d85', '#4e4466',
                 '#eccfb0', '#e06666', '#6fa8dc', '#0b5394', '#351c75',
                '#4b2328', '#a86800', '#395a57', '#9ab7c4', '#5e627b')

add_title <- function(x,title_str, subtitle_str = ''){
  x |>
  plotly::layout(title = list(text = paste0('<b>',title_str, '</b>','<br>','<sup>',subtitle_str,'</sup>')))
}


# Baixandos os dados da B3 (leva um tempinho, lembre-se de manter o cache para não deixar o TI da bolsa maluco)
# ajustesB3 <- futures_mget(first_date = '2006-01-01', last_date =  Sys.Date(), 
#                   do_cache = T)
# saveRDS(ajustesB3,'data/b3_ajustes.rds')

# Carregando os dados
ajustesB3 = readRDS('data/b3_ajustes.rds')

Obtidos os dados, vamos filtrar o data.frame para os contratos de DI1. Note que o fechamento vem em preço, será necessário transformar em taxa. Deixei a conta explicita.

di1_data = ajustesB3 |>
  filter(commodity == 'DI1') |>
  mutate(date_vencimento = rb3::maturity2date(maturity_code),
         date_vencimento_adj = bizdays::following(date_vencimento, 'Brazil/ANBIMA'), # Ajuste pois o vencimento pode cair em um feriado
         du_ate_vcto = bizdays(refdate, date_vencimento_adj, 'Brazil/ANBIMA'), # Calculando o número de dias úteis até o vencimento do contrato
         tx = ((100000/price)^(1/(du_ate_vcto/252))) -1) |> 
  select(refdate, date_vencimento, date_vencimento_adj, du_ate_vcto, price, tx) |>
   filter(du_ate_vcto > 0)

glimpse(di1_data)

Rows: 150,498
Columns: 6
$ refdate             <date> 2006-01-02, 2006-01-02, 2006-01-02, 2006-01-02, 2…
$ date_vencimento     <date> 2006-02-01, 2006-03-01, 2006-04-01, 2006-05-01, 2…
$ date_vencimento_adj <date> 2006-02-01, 2006-03-01, 2006-04-03, 2006-05-02, 2…
$ du_ate_vcto         <dbl> 22, 40, 63, 81, 124, 188, 249, 311, 373, 437, 499,…
$ price               <dbl> 98587.33, 97464.78, 96079.17, 95023.89, 92627.65, …
$ tx                  <dbl> 0.1769999, 0.1755996, 0.1734998, 0.1720998, 0.1684…

Interpolando as curvas e fixando as maturidades

Note que esta é uma estrutura dinâmica. O tempo passa, os contratos vencem. Para continuar o exercício, nós precisamos fixar as maturidades para trabalhar com taxas de 1,2,3 anos.

ultima_refdate_plot = di1_data |>
filter(refdate == max(refdate) | refdate %in% c('2023-09-01', '2023-01-03', '2022-10-04', '2020-04-01')) |>
mutate(refdate = as.factor(refdate)) |> 
ggplot() + 
   geom_point(aes(x = date_vencimento, y = tx, colour = refdate)) +
   geom_line(aes(x = date_vencimento, y = tx, colour = refdate)) +
   theme_minimal() +
   scale_y_continuous(labels = scales::percent) +
   scale_colour_manual('', values = paleta_cores) +
   labs(y = 'Taxa', x = 'Data de Vencimento do Contrato')

ggplotly(ultima_refdate_plot) # |>

  # add_title(title = 'Curva de Juros', subtitle_str = '% a.a.')

Para isso, é necessário interpolar as curvas — não se assuste, você está apenas ligando os vértices para cada data de referência. Assim, será possível pegar du_ate_vcto fixos de cada refdate. O pacote fixedincome facilita o processo. Vamos criar um objeto SpotRateCurve para cada data de referência, adicionar um método de interpolação e salvar apenas os vértices de interesse.

Quanto aos parâmetros: “discrete” é porque o método de capitalização utilizado é o discreto, “business/252” é porque nossa curva de juros é anualizada em 252 dias úteis e “Brazil/ANBIMA” para utilizar o calendário da ANBIMA para definir os dias úteis.

di1_data_l = di1_data |> 
  arrange(refdate, date_vencimento_adj) |>
  split(di1_data$refdate) |>
  purrr::map(function(x){
    curve = fixedincome::spotratecurve(x$tx, x$du_ate_vcto, refdate = unique(x$refdate),
                               'discrete', 'business/252', 'Brazil/ANBIMA')
    interpolation(curve) <- fixedincome::interp_naturalspline()
    curve
  })

curvaDI1 = tibble::tibble(refdate = names(di1_data_l), curvadi1 = di1_data_l) 

# Agora que temos um estrutura a termo interpolada para cada data de referência, conseguimos por exemplo, pegar a taxa de um ano para cada data de referência
dus = c(252, 252*2, 252*3, 252*4, 252*5, 252*6, 252*7, 252*8, 252*9, 252*10)

di1_constant_maturity = map_dfr(di1_data_l, 
          .f = function(x){
            tx = x[[dus]] |> as.numeric()
            tibble(refdate = x@refdate,
                   maturity = dus,
                   tx = tx)
          }
)

Para ilustrar o problema que acabamos de resolver, compare o gráfico abaixo com o primeiro gráfico que nós fixemos. Agora o gráfico não desloca mais, temos as taxas de 252 dias (1 ano), 504 dias (2 anos) e assim sucessivamente para cada ponto no tempo.

Com os dados organizados dessa forma conseguimos dar uma brincada e construir algumas visualizações interessantes.

# Gráfico com datas de referência selecionadas
p1 = di1_constant_maturity |> 
filter(refdate == max(refdate) | refdate %in% c('2023-09-01', '2023-01-03', '2022-01-03', '2021-01-04','2020-04-01', '2019-02-11')) |>
mutate(refdate = as.factor(refdate)) |>
    ggplot() +
    geom_point(aes(x = maturity, y = tx, colour = refdate), show.legend = F) +
   geom_line(aes(x = maturity, y = tx, colour = refdate), show.legend = F) +
   theme(legend.position = 'none') +
   theme_minimal() +
   scale_x_continuous(breaks = dus) +
   scale_y_continuous(labels = scales::percent) +
   scale_colour_manual('', values = paleta_cores) +
   labs(y = 'Taxa', x = 'DU')

ggplotly(p1) 
# Grafico interativo
p2 = di1_constant_maturity |>
    group_by(refdate_month = month(refdate), refdate_year = year(refdate)) |> # Pegando apenas o final de cada mês para não ficar muito pesado
    filter(refdate == max(refdate)) |>
    ungroup() |>
    # mutate(refdate = as.factor(format(refdate, '%b-%y'))) |>
    ggplot() +
      geom_point(aes(x = maturity, y = tx, frame = refdate), size = 2, colour = paleta_cores[1], show.legend = F)+
      geom_line(aes(x = maturity, y = tx, frame = refdate), size = 1, colour = paleta_cores[1], show.legend = F) +
      theme_minimal() +
      scale_x_continuous(breaks = dus) +
      scale_y_continuous(labels = scales::percent, n.breaks = 10) +
      # scale_colour_manual('', values = paleta_cores) +
      labs(y = 'Taxa', x = 'DU')

# Você pode criar um plotly com um seletor de datas
ggplotly(p2) |>
animation_opts(1000, easing = 'elastic', redraw = FALSE) |>
animation_button(x = 1, xanchor = 'right', y = 0, yanchor = 'bottom') 
#|>
# add_title(title = 'Curva de Juros', subtitle_str = '% a.a.')

# Todos os vértices ao longo do tempo
p3 = di1_constant_maturity |>
    mutate(maturity = as.factor(maturity)) |>
    ggplot() +
    geom_line(aes(x = refdate, y = tx, colour = maturity)) +
   theme_minimal() +
   scale_y_continuous(labels = scales::percent) +
   scale_colour_manual('', values = paleta_cores) +
   labs(y = 'Taxa', x = 'Data')

ggplotly(p3) 
#|>
# add_title(title = 'Curva de Juros', subtitle_str = '% a.a. por vértice')

# Superfície com o Plotly
di1_constant_maturity_xts = di1_constant_maturity |>
  pivot_wider(id_cols = refdate, names_from = maturity, values_from = tx)
di1_constant_maturity_xts = xts(di1_constant_maturity_xts |> select(-refdate), order.by = di1_constant_maturity_xts$refdate)

plot_ly(z = ~di1_constant_maturity_xts, y = index(di1_constant_maturity_xts), x = names(di1_constant_maturity_xts),
         colorbar = list(title = "Taxa")) %>%
  add_surface()  %>%
  layout(scene = list(
           legend = list('Taxa'),
           xaxis = list(title = "DU"),
           yaxis = list(title = "Data"),
           zaxis = list(title = "Taxa")
         )) #|>
  # add_title(title = 'Curva de Juros', subtitle_str = '% a.a., por maturidade e data de referência')

Análise de Componentes Principais (PCA)

Agora que temos tudo pronto, vamos seguir com a Análise de Componentes Principais. A ideia do PCA é que ele é capaz de para uma matriz de n x k, encontrar k-1 vetores independentes que contenham a maior parte da variabilidade encontrada nos vetores originais.

Complicando um pouco, o algoritmo busca pelos autovalores e autovetores da matriz de covariância ou correlação das séries, retornando k-1 vetores linearmente independentes.

# Vamos passar o df para wide e renomear os vértices para ficar mais claro
di1_constant_maturity_wide = di1_constant_maturity |>
  pivot_wider(id_cols = refdate, names_from = maturity, values_from = tx)
colnames(di1_constant_maturity_wide) <-  c('refdate',paste0(as.numeric(
  colnames(di1_constant_maturity_wide[-1])
  )/252, 'Y'))

psych::describe(di1_constant_maturity_wide[,-1]) |>
  knitr::kable(format = 'html',digits = 4)

	vars	n	mean	sd	median	trimmed	mad	min	max	range	skew	kurtosis	se
1Y	1	4315	0.1024	0.0329	0.1099	0.1046	0.0327	0.0219	0.1640	0.1421	-0.5429	-0.4565	5e-04
2Y	2	4315	0.1063	0.0292	0.1136	0.1079	0.0244	0.0323	0.1744	0.1421	-0.5063	-0.2413	4e-04
3Y	3	4315	0.1092	0.0264	0.1148	0.1104	0.0205	0.0421	0.1783	0.1362	-0.4246	-0.0904	4e-04
4Y	4	4315	0.1112	0.0245	0.1162	0.1122	0.0184	0.0497	0.1800	0.1303	-0.3712	-0.0121	4e-04
5Y	5	4315	0.1124	0.0232	0.1167	0.1133	0.0174	0.0542	0.1810	0.1268	-0.3291	0.0418	4e-04
6Y	6	4315	0.1134	0.0223	0.1173	0.1141	0.0173	0.0589	0.1814	0.1226	-0.2824	0.0782	3e-04
7Y	7	4315	0.1142	0.0216	0.1178	0.1148	0.0172	0.0625	0.1815	0.1190	-0.2382	0.1096	3e-04
8Y	8	4315	0.1148	0.0211	0.1181	0.1154	0.0169	0.0645	0.1815	0.1170	-0.2034	0.1640	3e-04
9Y	9	4315	0.1153	0.0206	0.1184	0.1157	0.0167	0.0657	0.1815	0.1158	-0.1689	0.2004	3e-04
10Y	10	4315	0.1156	0.0203	0.1185	0.1160	0.0167	0.0666	0.1815	0.1149	-0.1424	0.2451	3e-04

O gráfico a seguir mostra a correlação entre os vértices. Como esperado a partir da inspeção visual das séries, como mostra o gráfico a seguir a correlação é alta e positiva. É bastante provável que com poucos componentes principais consigamos expressar a maior parte da variabilidade das séries.

p = di1_constant_maturity_wide |>
    select(-refdate) |>
    corrr::correlate(quiet = T) |>
    corrr::rearrange() |>
    corrr::autoplot()
ggplotly(p)

Para quem já viu um pouco de finanças isso não é surpresa. Dentre as várias teorias que explicam a estrutura a termo, a Teoria das Expectativas diz que as taxas de juros de X anos nada mais é do que uma combinação da taxa de 1 ano e das taxas de 1 ano esperadas. Isso decorre de uma relação de não-arbitragem em que um investidor é indiferente entre investir em um título de 2 anos e investir em um título de 1 ano e depois reinvestir por mais 1 ano.

Para performar o PCA, temos a opção de tanto usar a matriz de covariância, como usar a matriz de correlação. Para essa aplicação especifica, vamos utilizar a matriz de correlação. Não queremos que o vértice de maior volatilidade domine sobre os demais. Esse problema fica bem mais evidente em outras aplicações em que as séries são bastante diferentes e com componentes idiossincráticos relevantes (e.g. taxas de câmbio).

# Vamos utilizar o pacote factoextra para performar o PCA
# Coloquei o processo numa função para extrairmos apenas o necessário
get_pca_results <- function(wide_data_frame){
  
  # PCA com normalização
  pca <- wide_data_frame |>
    select(-refdate) |>
    prcomp(scale. = T) # Scale = T utilizamos a matriz de correlação
  
  # Correlações
  pca_corr = factoextra::get_pca(pca)$cor |> as.data.frame.matrix() |>
    select(Dim.1:Dim.10) 
  
  # Contribuição de cada PC 
  pca_cotrib = pca |>
    broom::tidy(matrix = 'd') |>
    filter(PC <= 10) |>
    mutate(PC = as.character(PC)) 

  # Utilizando os loadings ou pesos para recuperar os componentes
  PCs = as.matrix(wide_data_frame |> select(-refdate)) %*% (pca$rotation*-1) |>
    as.data.frame.matrix() |>
    mutate(refdate = wide_data_frame$refdate) |>
    select(refdate, everything())
  
  list_results = list('prcomp.res'  = pca,
       'correlation' = pca_corr,
       'contrib' =  pca_cotrib,
       'PCs' = PCs)
  
  return(list_results)
}
di1_pca = get_pca_results(di1_constant_maturity_wide)

Analisando os componentes principais

Vamos aos resultados.

Primeiramente, sim, apenas 1 componente principal explica 97% da variabilidade em todos os 10 vértices que selecionamos da estrutura a termo. Os 3 primeiros componentes explicam praticamente 100%.

# Contribuição, Desvio Padrão por Componente principal
reactable(di1_pca$contrib |> filter(PC %in% as.character(c(1:5))),
          defaultColDef = colDef(format = colFormat(digits = 4)),compact = T, pagination = F, fullWidth = T) |>
  reactablefmtr::add_title('Contribuição e Desvio Padrão por Componente Principal', font_size = 16)

Contribuição e Desvio Padrão por Componente Principal

Mas o que tem em cada componente? Pois é, essa é a parte que mais me deixou curioso como um eterno aprendiz e iniciante em 999 coisas, embora para outras pessoas seja óbvio.

Para isso, vejamos as correlações dos componentes com as séries originais, bem como os loadings ou auto-vetores de cada componente. O loading nada mais é que a transformação linear que você aplica em cada série para obter o componente principal. Na prática, estamos construindo uma combinação linear de vértices. Cada combinação é linearmente independente da outra.

# Correlação entre o Componente Principal e o Vértice da Estrutura a Termo
reactable(di1_pca$correlation[,c(1:5)]*-1,
          defaultColDef = colDef(format = colFormat(digits = 2), style = reactablefmtr::color_scales(di1_pca$correlation*-1)), compact = T, pagination = F, fullWidth = T) |>
  reactablefmtr::add_title('Correlação entre os PCs e os Vértice do DI', font_size = 16)

Correlação entre os PCs e os Vértice do DI

O primeiro componente (PC1 ou Dim.1) tem alta correlação com todos os vértices e ela é bastante parecida em magnitude. De fato, como podemos ver no gráfico com os loadings, ele é uma combinação linear de pesos iguais de todos os vértices. Uma coisa que vou tentar mostrar posteriormente, é que, na prática, é como se estivessemos fazendo uma média de todos os vértices e extraindo um nível médio da curva naquela data de referência.

O segundo componente (PC2) tem correlação negativa com os vértices curtos e positiva com os vértices mais longos. É como se estivéssemos pegando a parte longa da curva e subtraindo da parte curta. Dito de outra forma, o PC2 é um fator que traz a inclinação da curva de juros.

Por fim, o PC3 ficou um pouco menos óbvio. O plot dos loadings tem quase um formato de U invertido, capturando um terceiro fator que é a curvatura da curva de juros.

# Loadings 
f=paste0(1:10, 'Y')
df_loadings = (di1_pca$prcomp.res$rotation*-1) |> as.data.frame()
df_loadings = df_loadings |> mutate(vertice = rownames(df_loadings)) |>
  select(vertice, everything()) |>
  mutate(vertice = factor(vertice, levels = f))

p=df_loadings |>
  pivot_longer(cols = -1, names_to = 'PC', values_to = 'loading') |>
  filter(PC %in% c('PC1', 'PC2', 'PC3')) |>
  ggplot() +
  geom_bar(stat = 'identity', aes(x = vertice, y = loading), fill = paleta_cores[6]) +
  # geom_hline(aes(yintercept = 0), size = .5, linetype = 'dashed') +
  theme_minimal() +
  facet_wrap(~PC) +
  labs(x = 'Vértice', y = '')
ggplotly(p) |>
  add_title(title = 'Loadings dos componentes principais', subtitle_str = '')

De fato, essas ideias estão bem documentadas na literatura de finanças. Com esses 3 fatores, devemos ser capazes de reproduzir toda a estrutura a termo da taxa de juros no tempo.

Comparando os componentes principais com média, inclinação e curvatura

Agora, nós vamos comparar: i) o primeiro componente principal (PC1) com uma métrica de nível da curva de juros; ii) o PC2 com uma medida de inclinação; e por fim, iii) o PC3 com uma métrica de curvatura.

Vamos definir a inclinação como sendo a difernça entre os juros de 10Y e os juros de 1Y. A curvatura será dada pelo dobro da média dos vértices de 4,5,6 e 7 anos subtraídos da soma dos dos juros de 1 e 10 anos. Já o nível é a média simples de todos os vértices.

#  Para isso vamos calcular os fatores nível, inclinação e curvatura com os dados originais
fatores_curva = di1_constant_maturity_wide |>
  mutate(nivel = rowMeans(di1_constant_maturity_wide[,-1]),
         inclinacao = `10Y`-`1Y`) |>
  rowwise() |>
  mutate(curvatura = 2*mean(c(`4Y`,`5Y`,`6Y`,`7Y`)) - (`1Y`+`10Y`)) |>
  ungroup() |>
  select(refdate, nivel, inclinacao, curvatura)

# Merge com os componentes principais
fatores_vs_pca = merge(di1_pca$PCs, fatores_curva, by = 'refdate')

fatores_vs_pca_long = fatores_vs_pca |>
  pivot_longer(cols = -1)

ay <- list(
  #tickfont = list(color = "red"),
  overlaying = "y",
  side = "right"
)

# PC1 vs. Nível (Média das Vértices)
plot_ly(data = fatores_vs_pca) %>%
  add_lines(x = ~refdate, y = ~nivel, name = "Nível",type = "scatter", mode = "lines",
            line = list(color = paleta_cores[5])) %>%
  add_lines(x = ~refdate, y = ~PC1, name = "PC1", yaxis = "y2",type = "scatter", mode = "lines", line = list(color = paleta_cores[4])) %>%
  layout(
    yaxis2 = ay,
    xaxis = list(title="Date", ticks=fatores_vs_pca$Date),
    yaxis = list(title = '')
  ) |> add_title(title_str = 'PC1 vs Nível da Curva de Juros')
# PC2 vs. Inclinação (10Y - 1Y)
plot_ly(data = fatores_vs_pca) %>%
  add_lines(x = ~refdate, y = ~inclinacao, name = "Inclinação",type = "scatter", mode = "lines", line = list(color = paleta_cores[5])) %>%
  add_lines(x = ~refdate, y = ~PC2, name = "PC2", yaxis = "y2",type = "scatter", mode = "lines",line = list(color = paleta_cores[4])) %>%
  layout(
    yaxis2 = ay,
    xaxis = list(title="Date", ticks=fatores_vs_pca$Date),
    yaxis = list(title = '')
  ) |> add_title(title_str = 'PC2 vs Inclinação da Curva de Juros',
                   subtitle = 'inclinação (eixo esq.), PC2 (eixo dir.)')
# PC3 vs. Curvatura (Média(2Y,3Y,4Y,5Y,6Y) - (1Y+10Y))
plot_ly(data = fatores_vs_pca) %>%
  add_lines(x = ~refdate, y = ~curvatura, name = "Curvatura",type = "scatter", mode = "lines",line = list(color = paleta_cores[5])) %>%
  add_lines(x = ~refdate, y = ~PC3, name = "PC3", yaxis = "y2",type = "scatter", mode = "lines",line = list(color = paleta_cores[4])) %>%
  layout(
    yaxis2 = ay,
    xaxis = list(title="Date", ticks=fatores_vs_pca$Date),
    yaxis = list(title = '')
  ) |> add_title(title_str = 'PC3 vs Curvatura da Curva de Juros',
                 subtitle = 'curvatura (eixo esq.), PC3 (eixo dir.)')

Como pode ser visto, os resultados são bastante satisfatórios para os dois primeiros componentes. No caso da curvatura, o fit já não é tão bom, sobretudo no passado. É possível melhorar fazendo 2*(3Y) - (1Y+10Y) — replicando um pouco mais fielmente os loadings obtidos para o PC3, entretanto imagino que a ideia de “curvatura” perderia um pouco o sentido. Sendo assim, vamos manter o resultado que encontramos. Um possível motivo para o shape dos loadings do terceiro componente não ter sido o esperado diante da literatura seja as discrepâncias de liquidez entre os vértices.

Enfim, essa é uma das intersecções mais interessantes que tive contato entre técnicas de Machine Learning e Finanças.

Enfim, obrigado se leu até aqui e se você encontrou algum erro ou alguma coisa muito estranha, pode me mandar mensagem ;).